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1 0 FREQUENCY ESTIMATION OF ELECTRO-ISLET-GRAPHY 



Field of the Invention 

The present invention relates generally to methods and systems for processing 
electro-physiological signals and, more particularly, concerns the detection and analysis of 
1 5 electrical activity related to blood glucose concentration. 

Background of the Invention 

Blood glucose level monitoring is of great importance for diabetics. 

Continuous monitoring of the glucose level can greatly reduce the medical complications, that 
20 are caused by metabolic imbalance. 

The Islets of Langerhans are located in the pancreas, and are responsible for 

the manufacture of insulin in the human body. An islet is a cluster of many cells. The Beta 

cells within the islets respond to glucose in bursts of electrical activity. The use of such islets, 

derived from the pancreas of donor animals, in a blood glucose monitor is disclosed in U.S. 
25 Patent 5,101,8 14. Electro-Islet-Graphy (EIG) is the measurement of the electrical activity of 

the islets of Langerhans. The present invention utilizes EIG to provide a continuous blood 

glucose level sensor. 

Studies demonstrate a clear correlation between the fundamental frequency of 

the EIG signal, and the glucose level in the medium surrounding the islet. Hence, the 
30 estimation of the fundamental frequency of Electro-Islet-Graphy is of significant practical 

value. 

The fundamental frequency is defined as the frequency of the "events" of the 
EIG. An "event" is believed to represent the synchronized electrical activity of the cells in the 
islet. By analogy to an ECG. an "event" in EIG is comparable to a the heart cycle (the PQRST 
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Although it is believed that EIG processing has not been performed by any 
entity but the owner the present patent application, one might attempt to detect the events 
directly, and then to calculate the fundamental frequency. The problem with this approach is 
5 that the shape and size of EIG "events" varies greatly, and therefore reliable and robust event 
detection is difficult to achieve. 

In accordance with the present invention, the fundamental frequency of EIG 
events is determined directly, without first detecting the individual events themselves. All 
these algorithms must use an analysis window containing more than one event. The present 
1 0 invention utilizes techniques similar to those utilized to detect and estimate pitch in speech 
processing 

The invention is based on two important biological discoveries. The first is 
that the EIG is generated by a functional pace maker. The second is that the EIG signal is 
quasi-periodic most of the time. Pitch detection algorithms are used, because of the 
1 5 essentially quasi-periodic nature of the EIG signal. By quasi-periodic we mean that (1) the 
intervals between successive events are not exactly identical, but may vary slightly and (2) the 
amplitude and shape of successive events may also exhibit some variance. 

Several pitch detection algorithms were tested. Three of them achieved good 
results: Autocorrelation, Segmented Autocorrelation and Harmonic Peaks analysis. The 
20 preferred embodiments of the invention focus on algorithms that are based on the 
Autocorrelation methods. 

The preferred version of the algorithm comprises the following steps: 



• Detection of the non-EIG signals. The signal may include non-EIG segments, such 
25 as artifacts and silences. The signal is scanned and the non-EIG segments are marked 

and ignored. 

• The signal is divided into overlapping analysis windows, each four 
seconds long and each has a 75% overlap with the adjacent windows. 
The analysis window contains more than one event. 
30 • A modified form of an Autocorrelation transform of the type used for 

pitch detection in speech processing is applied to a single analysis 
window. This step is repeated for each analysis window. 
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• The fundamental frequency is derived from the autocorrelation values 
of the analysis window. Usually the fundamental frequency is 
indicated by the largest autocorrelation value. This step is repeated for 
each analysis window. 

• A postprocessor is used to "smooth out" the results of all the individual 
analysis windows. The previously marked non-EIG segments ( 
artifacts and silences ) are added in this phase. 

To produce the modified autocorrelation transform, the existing autocorrelation 
based algorithm was adapted to EIG in the following ways: 

• An improved algorithm was devised for determining pitch from the 
autocorrelation values. The algorithm usually chooses the highest 
autocorrelation peak ( value ). In EIG we found that sometimes the 
true pitch is not represented by a peak, but rather by a valley between 
several adjacent peaks. An algorithm was devised to locate those 
cases, and to correctly estimate the pitch. We refer to this phenomena 
as a "volcano" shaped autocorrelation graph, because the center of the 
"mountain" is found on lower ground. 

• A voiced/unvoiced decision mechanism was adapted from speech 
processing. The "unvoiced" EIG segments were defined as a non- 
signals (undecided segments). A postprocessor was used to decide on 
the pitch of those undecided segments. Although unvoiced speech 
segments do exist, "unvoiced" EIG segments are a virtual non-signal, 
and do not really exist. 

• A special pre-processing algorithm was devised. The signal undergoes 
convolution so as to increase the width of the event. This is unlike 
speech pre-processing, which is aimed at enhancing the high amplitude 
portions and/or filtering out the formants. This preprocessing 
technique is referred to as "fattening". 

• A very long analysis window of 4 seconds was used, in order to find 
frequencies ranging from 0.25 Hz to 5 Hz. In speech it is customary 
to use a window of about 30 milliseconds, in order to find frequencies 
ranging from 80 Hz to 300 Hz. 
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The invention also contemplates the use of a Segmented Autocorrelation 
algorithm. This is considered to be the best method for EIG, but other methods are also 
adequate. Segmented autocorrelation is described in "Pitch detection of speech signals using 
Segmented Autocorrelation"/ 1. A. Atkinson, A. M. Kondoz, B.G. Evans/ Electronics Letters 
5 Vol. 31 No. 7 pp. 533-535/March 1995. 

Brief Description of the Drawings 

The foregoing brief description, as well as further objects, features, and 
advantages of the present invention will be understood more completely from the following 
1 0 detailed description of a presently preferred embodiment, with reference being had to the 
accompanying drawings, in which: 

FIG. 1 is a functional block diagram showing an overview of the process of the 

preferred embodiment; 

FIG. 2 is a flowchart describing the pre-detectors, the silence detector and the 

1 5 artifact detector; 

FIG. 3 is a flowchart describing fundamental frequency estimation of an 

individual analysis window; 

FIG. 4 is a flowchart expansion of one of the block 25 of FIG 3 and describes 
the decision algorithm of the fundamental frequency estimator; and 
20 FIG. 5 is a flowchart describing the post-processing phase algorithms, which 

combine information from neighboring analysis windows in order to correct local estimation 
errors. 

Detailed Description of the Preferred Embodiment 

2 5 The preferred embodiment of the invention will now be described in detail with 

reference to the accompanying drawings. First a schematic overview of the preferred 
embodiment of a process and system for fundamental frequency estimation of EIG is 
presented in FIG. 1. The algorithm comprises of the following steps: 

30 a. Detection of the non-EIG signals. The signal may include non-EIG segments, such as 
artifacts and silences. The signal is scanned and the non-EIG segments are marked and 
ignored. A silence detector 1 and an artifact detector 2 are used. An "artifact" is a 
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dominant interference noise signal, which severely corrupts the EIG signal. 



b. The signal is divided into overlapping analysis windows. Each window is 4 seconds 
long and has a 75% overlap with the adjacent windows. The analysis window must 

5 contain more than one event. It is contemplated that the analysis window could be as 

much as forty times the interval between successive events. 

c. The autocorrelation transform 4 is applied to a single analysis window. This step is 
repeated for each analysis window, pre-processing 3 is optional, and is not used in the 

1 0 preferred embodiment. 

d. The fundamental frequency is derived from the autocorrelation values of the analysis 
window 5. Usually the fundamental frequency is indicated by the largest 
autocorrelation value. A sureness measure 6 is computed for the estimated 

1 5 fundamental frequency of the analysis window. This step is repeated for each analysis 

window. 

e. A postprocessor 7 is used to "smooth out" the results of all the individual analysis 
windows. The previously marked non-EIG segments ( artifacts and silences ) are 

20 added in this phase. 

f . The output is the estimated fundamental frequency of the signal 8. 



The preferred embodiment of the invention includes software which preferably 
25 runs on a Pentium PC, using the Windows 95 operating system. The algorithm was 
implemented on the "Matlab" software by "The Mathworks Inc." The algorithm was written 
in the Matlab language, and runs within the Matlab program shell. It also requires the "Digital 
Signal Processing Toolbox" for Matlab, by "The Mathworks Inc." 

FIG 2. is a flowchart describing the pre-detectors. The silence detector ( 11 
30 , 12 ,13 ) involves two steps. In the first step 11 an amplitude related measure is computed for 
each second, and in the second step 12 the amplitude related measure is compared to a fixed 
threshold. "Sig" is the raw digital input signal. It is the raw signal recorded from the islet, 
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after an Analog to Digital conversion. In block 11 each one second window is examined. 
"Max(Sig)" is the maximum sample value within the examined one second window. 
"Min(Sig)" is the minimum of the samples 1 values within the examined one second window. 
As "Sig" is always a real number, the real maximum and minimum are computed. If the 
5 measure is below the threshold for more than 5 seconds, then the segment is classified as 
silence 13. The threshold depends on the measurement equipment and environment. It also 
depends on the digitizing scale. In the preferred embodiment the maximum amplitude of the 
EIG signal is approximately 4000 amplitude units, so a threshold of 100 amplitude units is 
used. An "amplitude unit" is the amplitude difference between two adjacent quantization 
1 0 levels of the Analog to Digital converter. In the preferred embodiment one "amplitude unit" 
corresponds to 0.25 Micro- Volt of the original electrical signal. "Original" means before 
amplification. 

The artifact detector ( 14, 15, 16, 17 ) involves two steps. In the first step 14 
the signal is compared to an adaptive threshold. The adaptive threshold is computed using a 
1 5 300 second long recording, by the following computation: 



A high reference is defined as the 99.833 percentile of the 
histogram of the amplitude values of the 300 second long 
recording. 



20 



A low reference is defined as the 0.167 percentile of the 
histogram of the amplitude values of the 300 second long 
recording. 



An amplitude reference is defined as the high reference minus 



25 



the low reference. 



A high threshold is defined as the high reference plus the 
amplitude reference. 



30 



A low threshold is defined as the low reference minus the 
amplitude reference. 
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The "histogram" used in block 14 is generated by simply sorting all the sample values 
("amplitudes") of a 300 second long recording. The sorting is done from the smallest value 
to the largest value. Based on those sorted values the two references are derived. The "high 
reference" is the 99.833 percentile of the sorted values. The "low reference" is the 0.167 
5 percentile of the sorted values. The term "histogram" refers to the sorting of the values, and 
can be replaced by the term "sorted amplitude values." 

The second step 15 involves checking whether there are samples in which the 
signal's amplitude is above the high threshold or below the low threshold. If such samples are 
found then they are classified as containing an artifact 16. The remaining samples 17 are not 
1 0 classified as an artifact or as a silence. 

FIG. 3 contains a flowchart describing the steps that are performed on each 
analysis window. Each analysis window is 4 seconds long. Overlapping analysis windows 
are used. There is a 3 second overlap between two successive windows. Analysis windows 
containing either silences or artifacts are not analyzed. The analysis windows are described 
1 5 in block 21, and are represented mathematically by the following notation: 
x(n) are the raw digital samples of the signal. We earlier 
referred to this raw input signal as "Sig." This is the same 
signal, only that here a more strict mathematical notation is 
used. The analysis window contains 4 seconds of signal. The 
20 sampling frequency ("FS ") in the preferred embodiment is 1 00 

Hz, therefore the analysis window contains 400 samples. This 
is marked by "N(Window length)=400." The first sample in 
the analysis window is x(l), where "7" is the index of the 
sample. Therefore the last sample in the analysis window is 
25 x(l+399)=x(l+N-l). Therefore the analysis window is 

described by the samples it contains from x(l) to x(l+N-l). 

Preprocessing 22 may be used, but it is not implemented in the preferred embodiment. 

A rectangular window 23 is applied to the samples, and a normalized biased 
30 autocorrelation transform 24 is computed. It is represented mathematically by the following 
equation: 
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q> (/n ) s * En-o m1 " 1 )w(n)][x(n+l+m)w(n+m)) $ 

where M w(n)'* is the windowing function." As "w(n)" is a rectangular window, as described 
5 in block 23, it is equivalent to "1" within the analysis window: 
w(n)= 1 0<;n<;N-l 
0 Otherwise 

The autocorrelation transform is further described in the book "Digital Processing of Speech 
Signals", L.R. Rabiner and R.W. Schafer, 1978, pp. 141-164. The fundamental frequency 

1 0 decision algorithm 25 is described in more detail in FIG. 4. 

We will now refer to FIG. 4 and later return to FIG. 3. 
FIG. 4 is a flowchart describing the fundamental frequency decision algorithm. 
Actually it decides on the fundamental period, from which the fundamental frequency can be 
derived. An improved algorithm was devised for determining the fundamental frequency from 

1 5 the autocorrelation values. In typical algorithms the highest autocorrelation peak (value) is 
usually chosen. In EIG, it was found that sometimes the true pitch is not represented by a 
peak, but rather by a valley between several adjacent peaks. An algorithm was devised to 
locate those cases, and to estimate the pitch correctly. We refer to this phenomenon as a 
"volcano" shaped autocorrelation graph, because the center of the "mountain" is found on 

20 lower ground. 

The input to the decision algorithm is a vector containing the autocorrelation 
coefficients 31. Steps 32 and 33 find the 5 highest extreme points within the allowed 
frequency range. A "volcano" effect detector 34 checks first whether one of the highest peaks 
is found in a suspected "volcano" area. If such a peak is found, it's index (/) is passed on in 

25 order to check the peak's amplitude 35 compared to the amplitude of the highest peak. If the 
result is positive then a "volcano" is detected 36 and the estimated fundamental period is half 
of the autocorrelation lag corresponding to the highest peak. The autocorrelation is 
represented by equation (1) above, m is the "autocorrelation lag." It is the lag between the 
two segments on which the correlation is computed. The first segment is x(n+l) to x(n+l+N- 

30 \m\-l) and the second segment is x(n+l+m) to (x+l+m+N-\m\-l). Also, see "Digital 
Processing of Speech Signals", L.R. Rabiner and R.W. Schafer, 1978, pp. 141-164. 

A pitch halving detection algorithm 37 first checks the ratio of the 



WO 99/60922 



PCT/US99/11973 



9 

autocorrelation coefficients' lag 37, and then their amplitudes 38. If one of the highest peaks 
fulfills both those conditions then it is chosen as the estimated fundamental period 39. If both 
the "volcano" detector and the pitch halving detector produced negative results, then the 
estimated fundamental period is assumed to be the lag of the highest autocorrelation peak 40. 
5 The output of the fundamental period decision algorithm of FIG. 4 is the estimated 
fundamental period. This output is returned to block 26 in FIG. 3. 

A "voiced" decision mechanism 26 (Fig. 3) is used to decide whether the 
fundamental frequency should be estimated in the current analysis window. If the maximum 
autocorrelation coefficient in the allowed lag range ( a lag of 20 to 333 samples ) is below 0.3, 

10 then the analysis window is classified as "unvoiced", and no fundamental frequency is 
estimated 27. In that case a sureness grade of zero is given to the analysis window. 

The voiced / unvoiced decision mechanism was adapted from speech 
processing. The "unvoiced" EIG segments were defined as non-signals ( undecided segments 
). The postprocessor will be used in later stages to decide on the fundamental frequency of 

1 5 those undecided segments. Note that unvoiced speech segments do exist, while "unvoiced" 
EIG segments are a virtual non-signal, and do not really exist. 

If the current analysis window is not classified as "unvoiced" then the 
estimation process continues, and a sureness grade is calculated for the current analysis 
window 28. A higher sureness grade indicates that there is a higher probability that the 

20 fundamental frequency was estimated correctly. For each analysis window several outputs 
29 are returned: 

1. The estimated fundamental period- this value is not returned if the analysis 
window is classified as "unvoiced"; 

2. A "voicedV'unvoiced" decision- specifies whether the fundamental period 
25 was estimated in the analysis window; and 

3. A sureness grade- this value is not returned if the analysis window is 
classified as "unvoiced". 

FIG. 5 describes the final stage of the estimation process: the postprocessor. 
The postprocessor combines all the acquired data 41 into a fundamental frequency estimate 
30 of the raw signal. The postprocessor operates on a 300 second long recording of the raw 
signal. The first step 42 is to convert all the estimated fundamental period values into 
estimated fundamental frequency ( "Pitch" ) values. The output of this stage is a fundamental 
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frequency estimate, in Hertz, for every "voiced" analysis window. The next step 43 inserts 
the previously detected silence segments 1 into the fundamental frequency estimation results. 
The fundamental frequency is assigned a value of 0 ( Hz ) in all windows containing only 
silences. All analysis windows containing an unlikely fundamental frequency estimate are 
5 marked as "unvoiced" in the next step 44. 

The next step 45 assigns a fundamental frequency ( "pitch" ) value to the 
"unvoiced" segments. The postprocessor loops once over all the analysis windows. It assigns 
Pitch values to all "unvoiced" windows ( windows containing artifacts are also treated as 
"unvoiced" ). Pitch values are derived from an average of the pitch values of the neighboring 

1 0 windows. Some minor data dependent modifications can be added in order to create a smarter 
averaging method, but those are minor changes which are very data dependent. A pitch 
halving / doubling error correcting mechanism is then applied 46. The algorithm loops once 
over all the analysis windows. If the Pitch of the current window is half or double the average 
of the neighboring windows, then the current Pitch is taken as a smart average of the 

1 5 neighbors. The term "smart averaging" refers to the average of an ensemble of numbers, after 
the extreme values are removed. For example: We can define a "smart average" of 8 numbers 
in the following way: 

1. Remove the maximal and minimal numbers and 

2. Compute the average of the remaining 6 numbers. 

20 This "smart averaging" method is used in order to achieve more stable and robust results. 

The final step 47 looks for an analysis window having a pitch which is 
significantly different than its neighbors. The algorithm loops once over all the analysis 
windows. If the Pitch of the current window, is 35% greater or smaller than the average of the 
neighboring windows, then the current Pitch is taken as a smart average of the neighbors. 

25 The final output 48 of the entire fundamental frequency estimation process is 

a fundamental frequency value for each analysis window. 

Although a preferred embodiment of the invention has been disclosed for 
illustrative purposes, those skilled in the art will appreciate that many additions, modifications 
and substitutions are possible without departing from the scope and spirit of the invention. 

30 For example, the principles arid applications of the invention are not limited to the particular 
biological phenomenon described. They could be utilized equally well, for example, for 
measuring an optical signal from living cells of a muscle or the electrical signals of neurons, 



WO 99/60922 POI7US99/1 1 973 

11 

or various signals produced by similar or other biological micro-structures. In particular, it 
is contemplated that the invention could find utility in any application where a cell or other 
biological micro-structure is moved outside its natural environment and used as a biological 
sensor. 



